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This  report  was  prepared  by  Don  W.  Kinsey  of  the  Aerodynamics  and 
Airframe  Branch,  Aeromechanics  Division,  Flight  Dynamics  Laboratory, 
Wright-Patterson  AFB,  Ohio.  The  work  was  performed  under  Work  Unit 
240410A1 . 

The  work  was  part  of  an  effort  to  develop  improved  computational 
capability  for  predicting  unsteady,  transonic  flows.  Modifications 
made  to  an  existing  computer  program  (LTRAN2)  are  described  and  some 
sample  results  are  provided. 
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I.  INTRODUCTION 


The  need  for  a  method  to  determine  the  unsteady,  transonic  flow 
conditions  over  an  airfoil  is  well  established.  The  possible  applica¬ 
tions  for  such  a  solution  method  include  buffet  and  flutter  studies 
for  fixed  wing  aircraft,  propeller  and  rotor  blade  studies,  turbine 
blade  and  guide  vane  designs  and  many  more.  The  equations  governing 
unsteady,  transonic  flow  of  a  viscous  fluid  (Navier-Stokes  equations) 
have  been  known  for  over  140  years.  However,  solving  the  equations 
for  any  practical,  time-accurate  transonic,  viscous  flow  over  ar¬ 
bitrary  airfoils  requires  vast  computer  resources.  Reducing  the 
problem  to  two  dimensions  still  leaves  a  challenging  task. 

Some  success  in  solving  attached  flow  problems  has  been  achieved 
by  making  simplifying  assumptions.  Flow  conditions  that  have  very 
thin  boundary  layers  can  be  closely  approximated  by  inviscid  equa¬ 
tions,  the  so-called  Euler  solutions.  Another  major  simplifying 
assumption  is  the  irrotational  flow  assumption.  This  allows  the 
definition  of  a  velocity  potential,  and  the  combination  of  the  equa¬ 
tions  of  motion  into  a  single  equation  referred  to  as  the  full 
potential  equation.  Finally,  the  full  potential  equation  can  be  ex¬ 
amined  and,  using  order-of -magnitude  arguments,  reduced  to  the  small 
disturbance  equation. 

Solution  of  the  two-dimensional,  transonic  small  disturbance 
(TSD)  equation  provides  a  relatively  quick  and  inexpensive  method  of 
predicting  unsteady  transonic  flow  over  thin  airfoils.  However,  flow 
conditions  which  produce  significant  boundary  layer  growth  and/or  sig¬ 
nificant  shock-boundary  layer  interactions  are  not  well  predicted  by 
the  inviscid,  irrotational  TSD  equations.  As  a  result,  calculations 
which  use  the  steady  small  disturbance  solution  as  a  starting  point 
for  unsteady  calculations  have  serious  inherent  errors. 

Information  about  the  steady-state  aerodynamics  on  the  surface  of 
an  airfoil  is  often  known  from  experimental  data  or  may  be  determined 
with  more  sophisticated  prediction  techniques.  If  a  new  geometry 


could  be  determined  so  that  the  TSD  solution  conforms  to  the 
prescribed  surface  aerodynamics,  a  much  more  satisfactory  steady-state 
starting  solution  would  be  established.  The  effects  of  boundary 
layers,  local  shock- induced  separation,  and  other  flow  conditions 
would  be  effectively  added  to  the  inviscid  solution.  Once  the  correct 
steady-state  conditions  are  produced  by  the  small  disturbance  proce¬ 
dure  applied  to  the  new  geometry,  subsequent  time  dependent  variations 
should  also  be  more  representative  of  the  real  flow  situation. 

A  numerical  procedure  for  determining  unsteady  aerodynamics  of 
two-dimensional,  thin  airfoils  using  the  TSD  equations  was  developed 

by  NASA/Ames  and  given  the  name  LTRAN21.  7  Several  versions  of  LTRAN2 
exist.  The  version  used  for  this  study  has  a  low-frequency  and  a 
high-frequency  option,  a  wind  tunnel  wall  option,  a  boundary  layer 
option  and  a  lag  entrainment  option.  An  implicit  algorithm  is  used  to 
solve  the  governing  equations  for  the  flow  about  thin,  two-dimensional 
bodies  in  oscillatory  motion,  plunge  motion,  or  an  oscillating  flap 
motion . 

The  normal  unsteady  solution  procedure  is  to  first  compute  the 
steady-state  potential  flow  about  the  body  at  a  given  initial  angle  of 
attack.  Then,  variations  in  the  potential  flow  about  the  initial  con¬ 
ditions  are  determined  as  a  function  of  time  for  the  type  of  motion 
def ined . 

The  objective  of  this  study  is  to  develop  a  procedure  that  allows 
the  user  to  obtain  the  initial  surface  potentials.  Once  the  poten¬ 
tials  are  defined  on  the  surface  of  the  airfoil,  a  modified  steady- 
state  solution  routine  is  employed  in  an  inverse  mode  to  obtain  the 
geometry  that  would  produce  the  target  surface  potentials.  This  new 
geometry  is  then  used  to  compute  unsteady  solutions  that  provide  much 
better  approximations  to  the  experimental  unsteady  characteristics. 
Experimentally  measured  surface  static  pressures  or  steady-state  solu¬ 
tions  from  higher  order  flow  solvers  are  both  excellent  sources  of 
surface  information  for  steady-state  flow  conditions. 

This  report  will  provide  a  description  of  the  governing  transonic 
small  disturbance  equations  and  the  boundary  conditions.  A  discussion 
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of  the  modifications  necessary  for  the  inverse  steady-state  procedure 
will  be  included  and  some  representative  results  will  be  shown. 


II.  GOVERNING  EQUATIONS  AND  SOLUTION  SCHEME 


Development  of  the  Small  Disturbance  Equations 

The  unsteady,  small  disturbance  f]ow  equation  is  derived  from  the 
unsteady  Euler  equations.  Assuming  irrotational  flow  (fluid  particles 
do  not  rotate  about  their  axis)  allows  the  Euler  equations  to  be  com¬ 
bined  into  a  single  equation  in  terms  of  a  velocity  potential,  4> .  For 
isentropic  flow,  relationships  exist  between  pressure  (p) ,  density  ( p ) 
and  speed  of  sound  (a)  so  that  p  and  p  do  not  appear  explicitly  in  the 
equation.  This  provides  the  full  potential  flow  equation: 


9^  d>  9  ?  2 

—  +ir(V$.V$)  *  V$.VF^)  =  a2  Tt 

at2  ot  z 


where, 


4  -  (7  -  1)[|?  ♦  (1/2)  (Vd>.V*  -  U2)] 


a^  =  free  stream  speed  of  sound. 


Further  simplification  of  the  potential  equation  is  possible  by 
restricting  the  analysis  to  cases  were  the  following  small  disturbance 
assumptions  are  valid; 


k,  1-M2,  and  r  ^  <<  1. 


Where,  k  is  the  reduced  frequency,  =  u  c/U^ 

c  is  the  chord  length  of  the  airfoil, 

U  is  the  undisturbed  freestream  velocity, 

00 

w  is  the  unsteady  oscillatory  frequency  of  the  airfoil, 
T  is  the  airfoil  thickness-to-chord  ratio, 


•V-V-'-.-VV- 
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M  is  the  freestream  Mach  number. 

w 

The  resulting  equation  can  be  written  in  terms  of  a  small  perturbation 
velocity  potential  (^)  as: 

Mtt  +  2B*xt  =  C,xx  +  ,yy.  (1) 

Where,  A  =  k2  M^  /r 

B  =  k  /r2/3 

c  =  [(lX  )/r2/3]  - (7+l)M^  0x. 

x  and  y  are  the  Cartesian  coordinates  in  the  freestream  and  vertical 
directions  respectively.  7  is  the  ratio  of  specific  heats.  m  is  nor¬ 
mally  2.0  (Spreiter  scaling)  but  may  vary  for  some  applications.  The 
subscripts  on  <p  represent  partial  derivatives  of  ^  with  respect  to  the 

2  /3 

subscripted  variable,  x  is  nondimensionalized  by  c,  y  by  c  r,  time 
by  w,'and  <f>  by  c  W^rJ^ 

LTRAN2  Description 

A  true  high  frequency  solution  would  require  the  solution  of 
equation  (1).  However,  LTRAN2  achieves  good  results  for  a  large  num- 

2 

ber  of  cases  by  dropping  the  <f>  term  of  0  [k  )  and  solving  the 

t/t/ 

equation3 


2  B  <p 


xt 


C  <f> 


XX 


yy 


(2) 


Both  the  low  frequency  and  the  high  frequency  options  in  LTRAN2  use 
equation  (2).  The  difference  between  the  low-frequency  and  high- 
frequency  calculations  of  LTRAN2  is  in  the  boundary  conditions  and  the 


expression  used  to  compute  pressure  coefficients.  The  high  frequency 
option  includes  time  dependent  terms  of  0[k]  <1  in  the  calculation  of 
the  pressure  coefficients  and  the  boundary  conditions  at  the  wake  and 

3  4 

downstream  boundary.  ’  A  summary  of  the  boundary  conditions  and 
pressure  coefficient  definitions  is: 


Variable 

Low  frequency 

High  frequency 

Pressure  coefficient 

C  =  -2  <fl  r2/3 
p  rx 

cp  = 

Wake  conditions 

«„]  =  0 

[flx^flt]  =  0 

Downstream  B.C. 

K  - 0 

flx^flt  =  0 

Upstream  B.C. 

fl  =  0 

II 

o 

Upper  and  lower  B.C. 

fl  =  0 
y 

o 

II 

where  the  brackets  []  are  used  to  signify  this  is  a  jump  condition 
across  the  wake . 

The  surface  tangency  boundary  condition  for  low  frequencies  is 
satisfied  by 

fly  =  fx(x>t)  for  y  =  0,  0<  x  >  1 

7  /?> 

where  y  =  T  f(x,t)  defines  the  airfoil  surface.  However,  both  the 
low-  and  high-frequency  options  of  the  version  of  LTRAN2  used  for  this 
study  solved  the  time  dependent  expression 


=  fX+  k  ft(X,t). 


(3) 


The  unsteady  part  of  LTRAN2  was  not  altered  for  the  results  in 
this  report.  Rather,  the  algorithm  for  determining  the  steady  solu¬ 
tion  was  modified  to  operate  in  both  an  inverse  or  a  direct  mode.  The 
modifications  necessary  for  the  inverse  procedure  are  discussed  in  the 
next  section.  But  first,  we  will  discuss  the  steady  equations  and  the 
steady  solution  algorithm. 

LTRAN2  Steady-State  Solver 

The  equations  of  interests  for  the  steady-state  case  are: 


C  <p  +  <p  —  0 
rxx  ryy 

(4) 

f  =  -0.5  C  r'7/i 

P 

(5) 

ip  =  f  (x)  . 

Ty  ) 

(6) 

Where  equations  (5)  and  (6)  are  defined  on  the  airfoil  surface. 
Equation  (6)  becomes  unbounded  at  the  leading  edge  of  airfoils  with  a 
finite  leading-edge  radius. 

The  finite  difference  solution  of  the  small  disturbance  equations 
are  computed  on  a  rectangular  grid  (Figure  1) .  The  x-axis  is  in  the 
freestream  direction  parallel  to  the  airfoil  chord  it  zero  angle  of 
attack  and  the  y-axis  is  vertical,  or  perpendicular  to  the  x  coor¬ 
dinate  . 

An  approximate  factorization  solution  algorithm  (AF2)  developed 

by  Ballhaus  et  al7.  is  used  in  LTRAN2  to  compute  the  steady-state  flow 
solution  for  equation  (4) .  Some  of  the  more  pertinent  aspects  of  that 
solution  scheme  are  given  here. 

A  finite  difference  approximation  for  equation  (4)  may  be  written 
by  defining  a  linear  operator  L  as: 


7 


.%  w  vvvxvvw*  w»  k.% .%  w  l-»  wv  _>  „>  „> 
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=  (C  <5  xx+  6  ^  -  b  =0. 


The  matrix  operator  L  contains  <5  and  <5  ,  which  are  second  central 

r  xx  yy’ 


difference  operators  in  x  and  y  respectively.  It  also  contains  C 
which,  for  purposes  of  analysis  is  considered  a  constant.  The  vector 
b  results  from  boundary  conditions  and  the  vector  ^  is  the  exact  solu¬ 
tion  to  the  finite  difference  matrix  equation  (7) . 

A  general  form  for  a  two-level  iterative  solution  procedure  for 
equation  (7)  is 


N  Dn  +  a  Rn  =  0 


where  Dn  =  -  (f>n  is  the  correction  to  (pn  after  n  iterations;  Rn 


L  <bn  =(C  6  *6  )in 

T  xx  yy 


b  is  the  residual  at  iteration  level  n;  and  a 


is  a  relaxation  parameter. 

The  proper  choice  of  N  is  the  key  to  convergence  acceleration. 

Two  guidelines  for  the  choice  of  N  have  been  established:  (1)  N  should 
be  chosen  to  resemble  L  as  closely  as  possible,  and  (2)  the  matrix 
operations  required  to  obtain  the  correction  vector  should  be  simple. 
An  expression  of  N  as  an  approximate  "factorization"  of  L  is, 


N  =  N  r  N 


where  factors  and  are  chosen  so  that  their  product  is  an  ap¬ 


proximation  to  L.  N]  and  N  must  also  be  chosen  so  that  the  solution 


procedure  is  stable.  Each  factor  may  have  a  bidiagonal  or  tridiagonal 


form. 


An  approximate  factorization  of  equation  (7)  that  has  the  neces¬ 


sary  characteristics  can  be  written  as: 


(P  -  C  $x)(/>  6x  -  5yy)Dn  =  p [  (C  5xx  ♦  5yy)*n  -  b] 


^  <y*V»r.  |T,  /-^AA.V*r.  -V  Sm  A/..*  /.A  A 

A  -V-.N  A  "AN  A  A  A  A  A  .NA  A  v  AA  A  A  A  AA  A  A  A  A  A  A  AAA 


m 


where  6  and  6  are  first  order  accurate  forward  and  backward  dif- 
x  x 

ference  approximations,  respectively;  and  p  is  a  parameter  to  be 
defined  later.  Note  that  equation  (10)  is  for  subsonic  flow.  For 
supersonic  flow  (C  >  0) ,  backward  differences  are  used  for  <5^  and  6 x 

Equation  (10)  can  be  rewritten  in  the  form 

Step  1:  (/J  -  C  £x)Fn  =  /? [ (C  5xx  *  5yy)^n  -  b] 


Step  2:  (p  6^  -  5yy)Dn  =  Fn  (U 

where  two  sweeps  through  the  grid  are  required  for  each  iteration. 

The  first  step  involves  a  matrix  solution  of  y  =  constant  grid  points 

and  defines  the  intermediate  function  Fn.  The  second  step  requires  a 
matrix  solution  for  each  x  =  constant  line  and  defines  the  new  poten¬ 
tials  at  each  grid  point. 


Acceleration  Parameters 


i 

In  the  approximate  factorization  equation  (10) ,  a  free  parameter 
P  was  introduced.  The  choice  of  /?  influences  the  convergence  rate 
and,  ideally,  a  new  value  of  P  that  would  provide  the  largest  reduc¬ 
tion  in  error  should  be  selected  after  each  iteration.  Unfortunately, 
no  procedure  for  defining  the  optimum  value  for  P  has  yet  been 
developed  for  transonic  flow  problems.  Therefore,  an  acceleration 
parameter  sequence  is  used  and  repeated  during  the  iteration  process. 
First,  estimates  are  made  of  the  values  of  p  required  to  minimize  the 
error  growth  for  the  high  ( p ^)  and  low  (/^)  frequency  extremes  of  the 


error.  Then,  P ^  and  p ^  are  incorporated  in  an  array, 


ft.  = 


t  = 


o 


ft  \ft  I  ft  1  (k  /(Q  1) 


i  n  i 


This  array  of  acceleration  parameters  is  used  sequentially  during  the 
course  of  iterating  the  solution.  An  eight-element  sequence  is  most 
often  used  (Q=8) .  Estimates  based  on  a  simple  linear  analysis 

provided  that  /J,  =  Ay  !  and  =  1 .  The  convergence  rate  may  be  im¬ 
proved  by  making  small  adjustments  to  the  values  of  and  / Jp  and 
from  the  use  of  other  types  of  sequences. 

Convergence  Criteria 

Some  criteria  are  needed  to  assess  the  degree  of  convergence 
after  each  iteration.  Typical  are  the  residual,  Rn,  the  correction, 
Dn,  the  aerodynamic  lift  coefficient,  C^,  and  the  number  of  supersonic 

points,  NSUP.  Values  for  Rn  and  Dn  approach  zero  as  the  solution  con¬ 
verges,  while  NSUP  and  C^  approach  a  fixed  value. 
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III.  INVERSE  PROCEDURE 


The  objective  of  the  direct  steady-state  solution  procedure  is  to 
compute  surface  pressure  information  from  known  surface  slope  informa¬ 
tion.  The  objective  of  the  inverse  method  is  to  estimate  surface 
slope  information  (and  thus  geometry)  from  known  surface  pressure  in¬ 
formation.  A  flow  chart  of  how  the  steady  solution  procedure 
currently  operates  is  shown  in  Figure  2. 

The  inverse  solution  procedure  required  several  changes  to  the 
direct  procedure.  The  major  difference  was  in  the  form  of  the  surface 
boundary  conditions.  In  the  direct  procedure,  the  surface  boundary 
conditions  are  Neumann  type,  where  the  normal  derivative  of  ^  was 
specified ; 


4>  =  f  (x) 

"  y  x v  ' 


where  f  is  the  local  slope  of  the  body.  In  the  inverse  procedure, 

the  boundary  conditions  on  the  surface  are  Dirichlet  type,  where  the 
value  of  the  variable  ^  itself  is  specified.  Both  of  these  boundary 
conditions  represent  a  well-posed  problem  that  has  a  unique  solution. 

Defining  Surface  Conditions 

The  first  step  in  the  inverse  procedure  was  to  define  the  0  dis¬ 
tribution  on  the  surface.  This  was  accomplished  by  converting  the 
input  pressure  coefficients  (known  from  experiment  or  higher  order 
computation)  into  surface  potentials  by  using  the  relationship  of 
equation  (5) 


K  =  <**)*  =  -°-5  T~? /3  c. 


where , 


(2  )f.  =  (*.  ^.)/(x.  x.)  =  -0.5  t~2/3  C 

v  x'rj  vyj  +  l  rj"  j  +  1  j  p. 


Solving  for  ^+1, 


1=  -0.5  t  C  (x.  .-x  )  +  j t- 
rj+i  Pj  j+i  r  r: 


j  is  the  grid  index  in  the  x  direction. 

Once  an  initial  estimate  for  <j>  at  x  =  0  was  determined,  the 
remaining  values  at  each  grid  point  (j)  along  the  body  were  computed 
from  equation  (13) .  The  results  presented  in  this  report  use  the 
value  of  ^  at  the  leading  edge  (x=0)  obtained  from  an  initial  direct 
solution  on  the  input  geometry.  Note  that  <f>  was  approximated  by  a 

forward  difference.  A  central  difference  approximation  for  0  was 

also  run,  but  there  was  no  obvious  improvement  in  the  results,  and  it 
had  the  disadvantage  of  requiring  two  initial  values  for  ^  at  the 
leading  edge. 

With  both  ^  and  ^  defined  on  the  surface,  it  was  a  simple  matter 
to  define  the  second  derivative  of  ^  in  the  x  direction  on  the  sur- 


0  =(0.  ,-2  <(> .  +  <f> .  .)/[0.5(x.  x.  -)]‘ 

rXXj  j-1  j  rj+l"  L  j+1  J-l 


Again,  several  different  approximations  for  are  possible,  but  the 

central  difference  appeared  to  work  best  for  the  cases  presented  in 
this  report. 

The  AF2  procedure  used  for  steady  solutions  in  LTRAN2  requires  an 
estimate  for  0  at  each  grid  point  on  the  surface  for  each  y  =  con¬ 
stant  sweep  iteration.  The  direct  method  computes  ^  from  the  known 

surface  slope  and  the  latest  approximation  for  <p  at  the  first  grid 
point  off  the  surface, 


where  I  is  the  surface  grid  line,  and  equation  (6)  is  the  relationship 

between  f  and  ^  .  Equation  (15)  is  essentially  a  forward  difference 
x  y 

approximation  for  the  second  derivative,  but  has  the  advantage  of 
needing  information  from  only  one  grid  level  away  from  the  surface  of 
the  body  (level  1+1) . 

The  inverse  solution  must  employ  another  procedure  for  estimating 
the  second  derivative  of  <j>  in  the  y  direction.  Second-  and  third- 
order  forward  difference  approximations  were  tried,  but  were 
unsatisfactory.  Very  large  gradients  can  exist  in  the  y  direction, 
particularly  near  the  leading  edge  of  the  airfoil,  and  the  forward 
difference  may  not  provide  a  good  approximation.  For  example,  the 
three  term  approximation  for  ^  ; 


I  -  2  *1  +  1  +  *I+2)/[°-5(yI  +  2~yI)]' 


has  errors  of  the  order 


^yI  +  l  *  1^  ^yyyj ^ 

which  may  not  be  a  small  value  when  large  gradients  exist.  Also,  the 
values  for  <p  on  the  surface  are  fixed,  while  the  values  away  from  the 
surface  are  evolving  during  each  iteration. 

8 

Following  the  example  of  Fung  and  Chung',  the  small  disturbance 
equation  itself  (equation  (4))  was  used  to  solve  for  ^  at  the  sur¬ 
face; 

<b  =  -C  <f> 

Tyy  rxx 


Substituting  for  C,  the  expression  becomes; 


tyy  ‘  -t,XldX  K  t,  I 


Note  that  this  approach  may  produce  inconsistent  results  just  above 
the  surface.  This  could  contribute  to  the  problems  with  convergence 
discussed  later  in  this  section. 


Equation  (16)  is  easily  solved  at  each  grid  point  on  the  surface 
during  the  y  =  constant  sweep  of  the  two  step  iteration  process. 

During  the  x  =  constant  sweeps,  the  code  was  modified  to  skip  the 
lines  of  code  that  redefine  the  values  of  ^  on  the  surface. 

Defining  the  New  Geometry 

After  the  inverse  process  has  converged,  the  potential  ^  is  known 
at  each  point  in  the  domain  defined  by  the  grid  system.  However,  the 
desired  information  is  the  shape  of  the  body  that  would  produce  this 
potential  distribution  when  input  to  LTRAN‘2  and  solved  in  the  direct 
mode.  Equation  (6)  provides  a  relationship  between  the  potentials  and 
the  surface  slopes; 

y'=  9y/9x  =  fx(x)  =  0y  =  (f[+1-f[)  /  (yi  +  r^l)  (17) 

Starting  with  y=0  at  x=0  (the  leading  edge  of  the  airfoil) ,  a 
Taylor  series  approximation  provides  an  estimate  for  y  at  the  next 
grid  point; 


2 

r.  =  y.  i  +  y‘.  (x.-x.  ,)  +  y''  (x.-x.  . )  /2  + 

J  J-l  J-l  J  J-l  J-l  J  J-l 


y-'t  (x.-x.  . )  /6  + . 
J-l  J  J-l 


(18) 


Where,  yj  ^  is  computed  from  equation  (17),  and 


with  similar  approximations  for  higher  derivatives. 

Geometries  were  computed  using  both  three  and  four  terms  of  the 
Taylor  series  approximation  and  with  forward,  centered  and  backward 
approximations  for  the  y' '  and  y’’”  terms.  The  best  results  were  ob¬ 
tained  using  the  first  three  terms  of  the  Taylor  series,  and  a  forward 
difference  approximation  for  y'. 


y  .  =  y  .  ,  +  <t  (x  x  .  . ) 
J  J-l  yj-l  J  J_1 


0. 5(x .  - 
J 


My  -  K  )] 

yj  yj-l  J  J 


or,  after  simplification, 


y.  =  y.  .  +  0.5(x.-  x.  ,)(^  +  <p  )  (19) 

'j  J-l  V  J  J"1  yj  yj~l 

A  new  geometry  is  created  by  successively  applying  equation  (19)  at 
each  grid  point,  starting  at  the  leading  edge  and  progressing  to  the 
trailing  edge.  Both  the  upper  and  lower  surfaces  are  constructed  in¬ 
dependently.  The  new  geometry  should  produce  a  pressure  distribution 
identical  to  the  target  pressures  when  input  to  the  direct  method  of 
LTRAN2.  When  this  new  geometry  and  the  resulting  steady-state  start¬ 
ing  solution  were  used  in  the  unsteady  mode,  a  much  better  comparison 
with  experimental  unsteady  data  was  obtained. 

Grid  Spacing  effects 

The  effect  of  grid  spacing  on  the  inverse  solution  was  inves¬ 
tigated.  Reduced  grid  spacing  in  the  x  direction  was  implemented  in 
the  leading  edge  region,  at  the  shock  location  and  at  the  trailing 
edge.  Both  the  leading  edge  and  the  shock  location  results  with  a 
finer  grid  spacing  were  worse  than  with  the  basic  LTRAN2  grid  dis¬ 
tribution.  However,  additional  grid  points  at  the  trailing  edge  did 
seem  to  help  smooth  the  pi essure  data  in  that  region.  Reducing  the 
grid  spacing  in  the  y  direction  at  the  body  had  a  significant  effect 
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on  the  results.  In  particular,  reconstruction  of  the  new  geometry, 
which  requires  an  accurate  estimate  of  ^  ,  was  much  better  with  the 

finer  grid  spacing.  Therefore  the  results  presented  in  this  report 
used  the  basic  LTRAN2  grid  in  the  x  direction  plus  two  additional 
points  at  the  trailing  edge  for  a  total  of  113  points,  50  of  which 
were  on  the  airfoil.  The  grid  extended  ±200  chord  lengths  with  Ax  = 
0.005  at  the  leading  edge  and  Ax  =  0.0028  at  the  trailing  edge.  The 
initial  grid  spacing  in  the  y  direction  was  0.001  for  the  inverse  pro¬ 
cedure  compared  to  the  normal  value  of  0.01  for  the  direct  solution. 
There  were  a  total  of  97  points  in  the  y  direction  which  extends  to 
±396  chord  lengths.  The  main  iteration  loop  used  1620  time  steps  for 
the  unsteady  calculations.  This  allowed  for  4.5  complete  cycles  (at 
360  steps  per  cycle)  so  that  any  transients  would  be  out  of  the  solu¬ 
tion  . 

Convergence  Characteristics 

Convergence  for  the  inverse  method  followed  much  the  same  pattern 
as  for  the  direct  method.  Thin  airfoils  in  subsonic  flow  converged 

rapidly  when  using  the  maximum  correction  (Dn) ,  or  change  in  lift 
coefficient  (C^)  as  the  convergence  criterion.  However,  the  maximum 

residual  (Rn)  would  decrease  only  about  three  orders  of  magnitude  for 
the  inverse  method.  Thick  airfoils  with  strong  shocks  had  maximum 
residual  reductions  of  only  three  orders  of  magnitude  for  both  the 
direct  and  the  inverse  methods.  The  maximum  correction,  number  of 
supersonic  points  (NSUP) ,  and  the  change  in  C^  all  indicated  good  con¬ 
vergence  after  about  500  iterations. 

Acceleration  Parameters  Used 

Many  different  combinations  of  acceleration  parameters  were 
tried.  The  number  of  parameters  in  the  sequence  (Q) ,  the  low  value 
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(/?^)and  the  high  value  (/J^)  may  be  changed  independently  for  each 
solution.  The  values  that  currently  exist  in  the  code  are: 


Process 

Direct  solution  on  input  geometry 
Inverse  solution 

Direct  solution  on  target  geometry 


There  was  little  attempt  to  optimize  these  parameters  for  the 
limited  cases  run  in  this  study.  Note  the  relatively  high  value  for 
in  the  inverse  method.  The  value  required  for  convergence  is  much 

larger  than  the  Ballhaus  et  al  estimate  of  1  for  the  direct  method, 

g 

but  is  consistent  with  the  findings  of  Chung  and  Fung. 

In  summary,  the  modifications  made  to  LTRAN2  for  this  report  in¬ 
clude  modifications  to  the  AF2  routine,  and  the  addition  of  four  small 
subroutines  to  accommodate  the  inverse  calculations.  Subroutine 
CPINPT  takes  the  target  surface  static  pressures  and  interpolates  to 
define  values  at  the  grid  point  locations.  SURPHI  computes  and  dis¬ 
tributes  values  of  ^  at  each  grid  point  on  both  the  upper  and  lower 
surface  (equation  (13)).  PYYC02  determines  the  second  derivative  of  <p 
in  the  y  direction  from  equation  (16) .  RECON  reconstructs  a  target 
geometry  based  on  the  results  of  the  converged  inverse  solution,  equa¬ 
tion  (19)  . 
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IV.  RESULTS  AND  DISCUSSION 


The  experimental  data  used  for  comparison  in  this  report  comes 

from  the  tests  of  Davis  and  Malcolm^  Two  airfoils  (NACA  64A010,  NLR 
7301)  were  tested  over  a  range  of  transonic  liach  numbers,  Reynolds 
numbers,  angles  of  attack  and  reduced  frequencies.  The  models  were 
oscillated  harmonically  about  50  percent  and  40  percent  chord  respec¬ 
tively.  Static  pressure  distributions  for  the  upper  and  lower 
surfaces  were  provided  for  each  test  condition  at  the  mean  angle  of 
attack  (aQ ) •  Dynamic  pressure  data  were  also  provided  at  several 

points  on  the  airfoil.  The  time-dependent  data  were  Fourier-analyzed 
up  to  the  fundamental  frequency  component;  where  the  fundamental  fre¬ 
quency  component  is  the  complex  number  that  indicates  the  magnitude, 
and  phase  shift  with  respect  to  the  input  motion.  The  real  part,  im¬ 
aginary  part,  magnitude  and  phase  of  the  complex  amplitudes, 
normalized  by  the  amplitude  of  the  input  motion,  were  tabulated  and 
are  used  for  comparison  with  LTRAN2  results  in  this  report.  Since 
magnitude  and  phase  are  related  to  the  real  and  imaginary  parts  by  the 
complex  relationship, 

C  +  i  C  =  C  e10, 
pR  PI  PM 

only  the  magnitude  and  phase  results  will  be  shown  in  this  report. 
64A010  Airfoil 

Table  1  lists  the  test  conditions  provided  in  this  report.  The 
first  test  condition  represents  the  most  severe  case  that  a  small  dis¬ 
turbance  procedure  could  reasonably  be  expected  to  handle.  The  64A010 
airfoil  is  10  percent  thick  and  has  a  moderately  rounded  leading  edge. 
A  significant  region  of  supersonic  flow  is  terminated  by  a  shock. 

The  airfoil  was  oscillated  at  a  reduced  frequency  of  0.202. 

Comparison  of  the  pressure  coefficients  from  the  experimental  test 
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(solid  line) ,  the  direct  solution  of  LTRAN2  on  the  original  64A010 
geometry  (dashed  line)  and  the  direct  solution  on  the  new  geometry 
that  was  determined  by  the  inverse  solution  (A  and  +) ,  is  presented  in 
Figure  3.  These  results  indicate  that  LTRAN2  predicts  the  steady  flow 
characteristics  for  this  case  with  reasonable  accuracy.  A  comparison 
of  the  original  geometry  in  Figure  4  (solid  line)  and  the  geometry 
determined  by  the  inverse  procedure  (+  and  A) ,  shows  very  little  dif¬ 
ference.  The  new  geometry  from  the  inverse  solution  has  some 
sharpening  of  the  leading  edge  and  som"  thickening  over  the  aft  part 
of  the  airfoil. 

The  time-dependent  results  of  LTRAN2,  for  both  the  original 
64A010  geometry  (BASIC  LTRAN2)  and  the  geometry  that  reproduces  the 
experimental  steady  results  (MODIFIED  LTRAN2) ,  and  the  experimental 
data  from  Ref.  9  are  compared  in  Figures  5  and  6.  Again,  there  ap¬ 
pears  to  be  little  difference  between  the  basic  results  and  the  new 
geometry  results,  and  both  compare  favorably  with  the  experimental 
results.  Some  representative  pressure  coefficient  plots  at  different 
stages  of  the  ocillatory  motion,  for  the  basic  airfoil  and  the  new 
geometry  respectively,  are  provided  in  Figure  7  and  8. 

NACA  0012  Airfoil 

The  second  test  case  was  the  NACA  0012  airfoil,  tested  at  Mach  = 

0.80,  and  1.86  degrees  angle  of  attack'.0  Reference  10  provides  no 
time-dependent  experimental  data,  but  the  0012  airfoil  is  the  ar¬ 
chetypal  test  airfoil  and  represents  a  reasonable  intermediate  shape 
between  the  64A010  and  the  NLR  7301  airfoils.  This  case  is  more 
severe  than  the  previous  case  because  the  NACA  0012  is  12  percent 
thick,  has  a  much  larger  leading  edge  radius  and  is  at  an  angle  of 
attack  of  1.86.  The  pressure  coefficient  comparisons  are  presented  in 
Figure  9.  The  direct  solution  for  the  NACA  0012  airfoil  is  obviously 
very  different  from  the  experimental  steady  data.  The  solution  for 
the  modified  geometry  agrees  quite  well  with  the  experimental  pressure 
data  except  for  a  small  deviation  at  the  shock. 


Theoretically,  the  inverse  procedure  should  produce  a  geometry 
that  would  exactly  reproduce  the  target  pressure  distribution.  In 
practice,  however,  it  has  been  found  that  pressure  data  in  the 
vicinity  of  a  strong  shock  is  very  erratic.  For  example,  the  symbols 
in  Figure  9  represent  the  pressure  results  after  1500  iterations  of 
the  inverse  procedure.  The  results  at  only  10  iterations  before 
(1490)  and  10  iterations  later  (1510)  show  that  an  almost  indistin¬ 
guishable  difference  in  geometry  produced  very  different  pressures 
(Figure  10) .  The  predicted  velocity  at  60  percent  chord  on  the  upper 
surface  varies  from  subsonic  to  supersonic  with  almost  imperceptible 
changes  in  geometry. 

The  geometry  provided  by  the  inverse  solution  using  the  Figure  9 
pressures  is  shown  in  Figure  11.  The  actual  NACA  0012  geometry  is 
much  thicker  and  has  a  much  more  rounded  leading  edge  than  the 
geometry  that  would  cause  LTRAN2  to  produce  the  target  pressure  dis¬ 
tribution.  The  "fishtail"  effect  on  the  aft  portion  of  the  new 
geometry  is  obviously  not  practicable  for  an  airfoil,  but  fortunately, 
the  small  disturbance  routine  only  requires  surface  slopes. 

The  unsteady  lift  and  pitching  moment  coefficients  from  the 
original  and  the  modified  geometries  are  compared  in  Figure  12.  The 
magnitude  and  phase  from  the  time-dependent  calculations  are  shown  in 
Figure  13.  A  look  at  the  predicted  pressures  at  various  points  during 
the  oscillations,  for  the  original  and  the  new  geometries  respec¬ 
tively,  is  given  in  Figures  14  and  15.  Several  cycles  were  computed 
to  allow  any  starting  transients  to  dissipate.  The  original  LTRAN2 
solution  is  virtually  invariant  with  angle  of  attack  and  produces  none 
of  the  phase  shift  and  complex  pressure  amplitudes  of  the  modified 
LTRAN2  procedure. 

NLR  7301  Airfoil 

Pressure  coefficients  for  the  NLR  7301  airfoil  are  shown  in 
Figure  16.  This  airfoil  is  16.5  percent  thick,  has  a  very  large  lead¬ 
ing  edge  radius  and  produces  strong  shocks  followed  by  separation  to 
the  trailing  edge  on  both  the  upper  and  lower  surfaces.  As  with  the 


NACA  0012  example,  the  direct  solution  for  the  basic  NLR  7301  geometry 
is  completely  wrong.  The  solution  for  the  modified  geometry  is  very 
good  over  the  first  50  percent  and  the  last  35  percent  of  the  airfoil. 
However,  the  small  disturbance  procedure  predicts  too  severe  a 
velocity  decrease  across  the  shock.  Attempts  to  improve  the  situation 
by  changing  grid  spacing  proved  to  be  fruitless.  It  may  well  be  an 
insurmountable  characteristic  of  the  small  disturbance  approximation. 

Indeed,  Steger  and  Klineberg1 1  show  similar  effects  in  their  attempts 
to  use  the  small  disturbance  equations  in  the  design  mode. 

The  original  and  inverse  airfoil  geometries  are  shown  in  Figure 
17.  Again,  the  large  leading  edge  radius  of  the  basic  airfoil  leads 
to  a  fishtail  effect  for  the  new  geometry. 

The  direct  solution  on  the  original  geometry  results  are  shown  in 
Figures  18  and  19.  The  original  LTRAN2  results  show  very  little 
change  in  the  time-dependent  pressure  profiles  and,  therefore,  very 
little  change  in  the  complex  pressure,  magnitude  and  phase.  The 
Figures  18  and  19  results  are  for  k  -  0.20.  A  k  =  0.05  case  was  also 
computed  but  the  results  were  as  bad  as  or  worse  than  the  k  =  0.20 
results . 

The  k  =  0.20  results  for  the  new  geometry  from  the  inverse  proce¬ 
dure  are  presented  in  Figures  20  and  21.  The  LTRAN2  results  compare 
favorably  with  the  experimental  results,  though  the  magnitude  of  the 
peak  values  are  different.  The  reason  for  the  differences  is  not  com¬ 
pletely  understood  but  is  very  likely  associated  with  the  shock 
strength  errors. 

Results  for  the  k  =  0.05  case  are  presented  in  Figures  22  and  23. 
As  with  the  k  =  0.20  data,  the  form  of  the  results  compares  with  the 
form  of  the  experimental  results,  however,  the  agreement  is  not  quite 
as  good  as  for  k  =  0.20.  Pressure  coefficient  plots  in  Figure  23  in¬ 
dicate  that  the  solution  in  the  neighborhood  of  the  shock  strongly 
influences  the  unsteady  pressure  results.  Additional  results  and 
recommendations  for  further  work  are  presented  in  Section  V. 
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V.  CONCLUSIONS 


The  steady-state  solution  subroutine  in  the  unsteady  small  dis¬ 
turbance  solver  LTRAN2  has  been  modified  to  operate  in  an  inverse  or 
"design"  mode.  "Design"  is  written  with  quotation  marks  because  there 
has  been  no  attempt  to  produce  a  realistic  airfoil  shape.  Rather,  the 
objective  of  this  inverse  procedure  is  to  provide  a  geometry  (actually 
a  set  of  upper  and  lower  surface  slopes)  that  causes  LTRAN2  to  match 
the  known  steady-state  surface  pressure.  The  unsteady  results  from 
the  new  geometry  should  be  a  much  better  approximation  to  the  ex¬ 
perimental  unsteady  results.  The  purpose  is  to  expand  the  useful 
range  of  the  small  disturbance  solver  to  both  airfoils  and  flow  condi¬ 
tions  that  violate  the  inviscid,  irrotational ,  small  disturbance 
assumptions . 

Three  separate  airfoils  were  evaluated  with  this  new  procedure. 
First,  an  NACA  64A010  airfoil,  tested  at  Mach  =0.80  and  near  zero 
angle  of  attack  was  evaluated.  This  condition  is  representative  of 
the  most  severe  case  that  a  small  disturbance  procedure  could  be  ex¬ 
pected  to  handle.  Indeed,  there  was  little  difference  between  the 
r  s  obtained  from  LTRAN2  run  in  the  normal  mode  and  the  results 
ol  ned  from  the  modified  geometry  that  produced  steady-state  pres¬ 
sure  distributions  that  exactly  matched  experimental  data.  The 
unsteady  pressure  results  from  both  approaches  agreed  well  with  ex¬ 
perimental  unsteady  results. 

Next,  an  NACA  0012  airfoil  at  Mach  =  0.80  and  a  =  1.86  was  inves¬ 
tigated.  The  steady-state  flow  results  indicated  clearly  that  the 
LTRAN2  solution  was  not  representative  of  the  experimental  data. 
Implementing  the  inverse  procedure  provided  a  geometry  that  did  repre¬ 
sent  the  experimental  results  when  solved  by  LTRAN2.  A  problem  of 
extreme  sensitivity  of  the  pressure  coefficients  at  the  shock  was 
revealed.  LTRAN2  over  estimates  the  velocity  reduction  at  the  foot  of 
strong  shocks.  There  were  no  experimental  unsteady  results  for  com¬ 
parison,  but  the  results  from  the  modified  geometry  were  much  more 
reasonable  that  the  results  from  the  basic  NACA  0012  geometry. 


The  most  severe  test  case  was  an  NLR  7301  airfoil  at  Mach  -0.80 
and  near  zero  angle  of  attack.  This  16.5  percent  thick  supercritical 
airfoil  had  a  very  large  leading  edge  radius.  The  flow  contained 
strong  shocks  on  both  surfaces.  Again,  the  steady-state  solution  on 
the  basic  geometry  was  much  different  from  the  experimental  steady- 
state  results.  The  inverse  procedure  produced  a  geometry  that 
resulted  in  good  agreement  with  the  experimental  steady-state  results 
except  for  the  "overshoot"  problem  at  the  foot  of  the  shock. 

The  unsteady  results  are  compared  for  two  different  reduced  fre¬ 
quencies;  k  =  0.20  and  k  =  0.05.  The  trends  of  the  computed  k  =  0.20 
results  compared  well  with  the  trends  of  the  experimental  results. 
However,  the  magnitudes  of  some  peak  values  were  significantly  dif¬ 
ferent.  The  k  =  0.05  results  did  not  compare  quite  as  well  with  the 
experimental  data  as  the  k  =  0.20  results. 

The  k  =  0.05  results  from  the  NLR  7301  airfoil  highlight  two 
problem  areas  for  the  inverse  procedure.  Probably  the  most  serious 
problem  results  from  the  fact  that  the  experimental  shock  location 
moves  fore  and  aft  on  the  airfoil  as  the  airfoil  oscillates.  However, 
the  modified  geometry  has  surface  slopes  that  induce  the  shock  for  the 
steady-state  solution,  and  this  geometry  remains  fixed  as  the  airfoil 
oscillates.  As  a  result,  the  geometry  is  trying  to  force  a  rapid 
change  in  surface  velocity  (and  thus  pressure)  at  the  wrong  (except, 
perhaps  twice  during  each  oscillation)  chord  location. 

The  second  problem  area  has  been  mentioned  before  and  is  the  ten¬ 
dency  for  LTRAN2  to  overshoot  the  velocity  drop  across  the  shock.  The 
experimental  data,  at  least  for  the  NLR  7301  case  investigated  here, 
indicate  a  much  less  severe  velocity  drop  across  the  shock  than  LTRAN2 
predicts.  The  overshoot  appears  to  be  an  unavoidable  result  of  the 
LTRAN2  procedure.  Some  minor  modifications  to  the  code  were  tried  but 
they  did  not  improve  the  overshoot  problem. 

Before  using  this  procedure  on  a  routine  basis,  a  parametric 
study  of  the  acceleration  parameters  should  be  completed.  It  would 
probably  be  very  useful  to  develop  a  procedure  to  vary  the  accelera¬ 
tion  parameters  as  the  solution  evolves.  The  solutions  presented  in 
this  report  were  obtained  on  a  CRAY  X-MP  12,  with  run  times  of  about  4 
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minutes,  and  required  about  156,000  words  of  memory.  However,  there 
was  no  attempt  to  improve  the  vectorization  of  the  code  or  to  decrease 
the  memory  requirements.  Significant  reductions  in  both  run  time  and 
memory  requirements  should  be  possible.  A  more  difficult  but  much 
more  rewarding  improvement  would  result  from  solving  the  overshoot 
problem  at  the  shock.  Last,  but  certainly  not  least,  the  development 
of  a  practical  way  of  allowing  the  geometry  modif ications  to  move  with 
the  shock  during  the  unsteady  oscillations  should  reap  tremendous 
benefits . 
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Table  1 


Test  Conditions  for  LTRAN2  Solutions 


Airfoil 

Mach 

k 

Xo 

R 

n 

% 

a 

CL 

NACA  64A010 

0.796 

0.202 

0.50 

1.3xl07 

-0.21 

1.09 

-0.029 

NACA  0012 

0.80 

0.20 

0.50 

9.0x10® 

1.86 

1.0 

0.299 

NLR  7301 

0.807 

0.20 

0.398 

1.2xl07 

0.36 

0.49 

0.047 

NLR  7301 

0.807 

0.05 

0.404 

1.2xl07 

0.36 

0.49 

0.047 

W 
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PRESSURE  CDEPriClENTS 


▼1  LOWER  SURFACE 


•UPPER  SURFACE 


-  EXPERIMENTAL  PRESSURE  COEFFICIENTS 

-  LTRAN2  SOLUTION 

+  A  INVERSE  LTRAN2  SOLUTION 


Figure  3  Pressure  Coefficient  Comparison,  NACA  64A010,  M=0.796 

a  =  -0.21,  R  =  1.3xl07 
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PRE5SURE  COCrriCIENTS 


-  LTRAN2  SOLUTION 

A  INVERSE  LTRAN2  SOLUTION 

—  EXPERIMENTAL  RESULTS 


Cl  vs 


Basic  LTRAft 2 


Modified  LTRAN2 


Basic  LTRAN2 


-  EXPERIMENTAL  PRESSURE  COEFFICIENTS 

-  LTRAN2  SOLUTION 

A  INVERSE  LTRAN2  SOLUTION 
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LIST  OF  SYMBOLS 


C^  or  CL 

C  or  CM 

m 

C  or  CP 
P 


2  2  2  /3 

Small  disturbance  equation  constant  =k  Mw  /r 
Speed  of  sound. 

Approximate  factorization  scheme  2. 

Angle  of  attack. 

2  2/3 

Small  disturbance  equation  constant  =k  Mb  /r 

Boundary  condition  vector  in  AF2. 

TSD  equation  = [ (l-M^)/x2 /3 ] - (7+l)l£  ^ 

Airfoil  chord  length. 

Lift  coefficient. 

Pitching  moment  coefficient  about  0.25c. 
Pressure  coefficient. 


Pressure  coefficient  corresponding  to  sonic  flow. 
Real  part  of  complex  pressure  coefficient. 


Imaginary  part  of  complex  pressure  coefficient. 


LTRAN2 


2  2l/2 

Magnitude  of  complex  pressure  coefficient  =  (C^  +  C^)  7 

Function  that  defines  the  airfoil  geometry. 

Intermediate  result  for  the  two  step  iteration  process. 
Reduced  frequency  =w  c/U^ 

Linear  operator. 

Unsteady  small  disturbance  solution  procedure. 

Mach  number. 

Exponent  for  Mach  number  (real  number  between  1  and  2) . 
Operator  used  in  the  AF2  development. 

Number  of  supersonic  points. 


LIST  OF  SYMBOLS 


Of  the  order  of  what  follows  in  brackets 

Number  of  parameters  in  the  acceleration  sequence. 

Residual  in  the  AF2  procedure. 

Reynolds  number  (per  foot) . 

Time  variable. 

Freestream  velocity. 

Cartesian  coordinate  in  the  streamwise  direction. 

Chord  location  about  which  the  forced  oscillations  occur 

Chord  location  in  fraction  of  chord. 

Cartesian  coordinate  perpendicular  to  the  freestream. 
Vertical  coordinate  in  fraction  of  chord. 


Angle  of  attack. 

Acceleration  parameter  used  in  AF2  solution  algorithm. 
Angle  of  attack  about  which  the  oscillations  occur. 

Multiplier  for  the  trig  function  describing  the 

oscillations . 

Upper  bound  for  the  acceleration  parameters. 

Lower  bound  for  the  acceleration  parameters. 

Thickness  to  chord  ratio  for  the  airfoil. 

Partial  derivative  operator  (with  a  spatial  or  temporal 
subscript) . 

Oscillatory  frequency. 

Velocity  potential . 

Small  perturbation  velocity  potential. 
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LIST  OF  SYMBOLS 


7  Ratio  of  specific  heats. 

a  Relaxation  parameter  for  the  AF2  solution  algorithm. 

9  Phase  angle  for  complex  pressure  coefficient  =  Tan  '(C  /C  ) 

PR  PI 


Subscripts 

«  Infinity  or  far  field  conditions, 

x  Derivative  with  respect  to  the  x  coordinate, 

y  Derivative  with  respect  to  the  y  coordinate, 

t  Derivative  with  respect  to  time. 

i  Index  of  the  computational  grid  in  the  y  direction. 

I  i  index  at  the  airfoil. 

j  Index  of  the  computational  grid  in  the  x  direction. 

J  j  index  at  the  airfoil  leading  edge, 

k  Index  for  the  acceleration  parameters. 


Superscripts 

n  Iteration  counter  for  the  AF2  solution  procedure. 

•*  Forward  difference. 

«■  Backward  difference. 

Derivative  with  respect  to  the  x  coordinate. 
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